Use marmap::getNOAA.bathy() to load in New Zealand bathymetry data directly from NOAA.

nz <- getNOAA.bathy(
  lon1 = 162,
  lon2 = 180,
  lat1 = -33,
  lat2 = -50,
  resolution = 1
) 
## Querying NOAA database ...
## This may take seconds to minutes, depending on grid size
## Building bathy matrix ...

Use sf::read_sf() to read the shapefile containing the New Zealand coastline.

nz_coast <- sf::read_sf(
  here(
    "data",
    "nz-coastlines-and-islands-polygons",
    "nz-coastlines-and-islands-polygons-topo-150k.shp"
    )
  )

Base R plot

blues <- c("lightsteelblue4", "lightsteelblue3",
           "lightsteelblue2", "lightsteelblue1")

greys <- c(grey(0.6), grey(0.93), grey(0.99))

plot(nz, image = TRUE, land = TRUE, lwd = 0.03,
     bpal = list(c(0, max(nz), greys),
                 c(min(nz), 0, blues)))
# Add coastline
plot(nz, n = 1, lwd = 0.4, add = TRUE)

Using the marmap::autoplot.bathy() function to plot with ggplot.

autoplot.bathy(nz, geom=c("tile","contour")) +
    scale_fill_gradient2(low="dodgerblue4", mid="gainsboro", high="darkgreen") +
    labs(y = "Latitude", x = "Longitude", fill = "Elevation") +
    coord_cartesian(expand = 0)+
    ggtitle("A marmap map with ggplot2")
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

Convert raster object to dataframe for ggplot.

nz %>% 
as.raster() %>% 
  rasterToPoints() %>% 
  as.data.frame() ->
  nz_df

Plot New Zealand Coastline with bathymetric and elevation data.

ggplot() +
  geom_raster(data = nz_df 
              # ## Uncomment below to only get bathymetry
              # %>% 
                 # filter(layer <= 0)
              , 
              aes(x = x, y = y, fill = layer)) +
  # ## I think it looks better without the contours on land, if any at all
  # ## Uncomment first and last two lines below to get all contours
  # ## Uncomment all below lines to get only seafloor contours 
  # geom_contour(data = nz_df
  #              %>%
  #                filter(layer <= 0)
  #              ,
  #              aes(x = x, y = y, z = layer), color = "grey20", size = .1) +
  scale_fill_viridis_c(option = "A") +
  geom_sf(data = nz_coast, fill = NA, color = "black", size = .15) +
  coord_sf(xlim = c(162, 180), ylim = c(-33, -50), expand = c(0, 0)) +
  labs(fill = "Elevation\n&\nBathymetry",
       x = "Longitude",
       y = "Latitude", 
       title = "New Zealand",
       subtitle = "Digital elevation (1 arc minute resolution)") +
  theme_classic()

Convert from bathy object to RasterLayer Object

nz %>% 
  as.raster() ->
  nz_raster 

Plot Raster onto interactive leaflet map

leaflet() %>%
  addProviderTiles(providers$Esri.OceanBasemap, group = "Ocean Base Map") %>% 
  addRasterImage(nz_raster, group = "Bathymetry") %>% 
  addPolygons(data = nz_coast, weight = 1, opacity = 1, fill = FALSE) %>% 
  addLayersControl(
    baseGroups = c("Ocean Base Map"),
    overlayGroups = c("Bathymetry"),
    options = layersControlOptions(collapsed = FALSE)
  )
nz_big <- getNOAA.bathy(
  lon1 = -180,
  lon2 = 180,
  lat1 = 90,
  lat2 = -90,
  resolution = 7
) %>% 
  as.raster()
## Querying NOAA database ...
## This may take seconds to minutes, depending on grid size
## Building bathy matrix ...
leaflet() %>%
  addProviderTiles(providers$Esri.OceanBasemap, group = "Ocean Base Map") %>% 
  addRasterImage(nz_big, group = "Bathymetry") %>% 
  addLayersControl(
    baseGroups = c("Ocean Base Map"),
    overlayGroups = c("Bathymetry"),
    options = layersControlOptions(collapsed = FALSE)
  ) %>% 
  setView(lng = 170, lat = -40, zoom = 4)